/*
    KMC Simulation
    Copyright (c) 2012 Filip Rozbořil

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "GridField.h"
#include "Grid.h"
#include "Utils.h"
#include "KmcSettings.h"
#include "Exception.h"
#include "constants.h"

bool operator==(const GridField& gf1, const GridField& gf2)
{
	//Trace("[OPERATOR==] Determining equality of " << gf1.mX << "," << gf1.mY << " and "
	//	<< gf2.mX << "," << gf2.mY);
	return gf1.mX == gf2.mX && gf1.mY == gf2.mY;
}

GridObject GridObject :: None;

const double Grid :: mBoltzmannConst = 8.617e-5; // [eV/K]

/*
	Constructs a Grid object. Extracts grid parameters and energies from the settings object.
*/
Grid :: Grid(KmcSettings& sett) : mEnergies(Site::Count), mLastHash(0)
{
	sett.getValue("Rows", mRows);
	assert(mRows > 4 && Utils::isPowerOfTwo(mRows));
	
	sett.getValue("Cols", mCols);
	assert(mCols > 4 && Utils::isPowerOfTwo(mCols));
	
	mElements = (uint64_t)mRows * (uint64_t)mCols;
	
	sett.getValue("Coverage", mCoverage);
	assert(mCoverage > 0.0 && mCoverage <= 0.5);
	
 	double totalAtoms = mCoverage * mElements + 100;
 	
 	sett.getValue("Temperature", mTemperature);
	assert(mTemperature > 0.0);
 	
 	mObjectReservation = (uint32_t) (1.1 * totalAtoms);
 	mVectorReservation = (uint32_t) (0.15 * totalAtoms);
 	
 	mAtoms.reserve(mObjectReservation);
 	mObjects.rehash(mObjectReservation);
 	
 	Trace("Initializing energies");
 	
 	sett.getValue("BarrierVertical", mBarrierVertical);
 	assert(mBarrierVertical > 0.0);
 	
 	sett.getValue("BarrierHorizontal", mBarrierHorizontal);
 	assert(mBarrierHorizontal > 0.0);
 	
 	sett.getValue("LogFrequencyVertical", mFrequencyVertical);
 	mFrequencyVertical = pow(10.0, mFrequencyVertical);
 	
 	sett.getValue("LogFrequencyHorizontal", mFrequencyHorizontal);
 	mFrequencyHorizontal = pow(10.0, mFrequencyHorizontal);
 	
 	double energyKink, energySingle, energyDimer, energyLateral;
 	
 	sett.getValue("EnergyKink", energyKink);
 	sett.getValue("EnergySingle", energySingle);
 	sett.getValue("EnergyDimer", energyDimer);
 	sett.getValue("EnergyLateral", energyLateral);
 	
 	mEnergies[Site::Free] = 0.0;
 	mEnergies[Site::Kink] = energyKink;
 	mEnergies[Site::Single] = energySingle;
 	mEnergies[Site::DoubleKink] = 2 * energyKink;
 	mEnergies[Site::Lateral] = energyLateral;
 	mEnergies[Site::Dimer] = energyDimer;
 	mEnergies[Site::DimerSingle] = energyDimer + energySingle;
	mEnergies[Site::DimerKink] = energyDimer + energyKink;
	mEnergies[Site::DimerDoubleKink] = energyDimer + 2 * energyKink;
}

/*
	Calculates the frequency of hopping from site "from" to site "to" in direction "dir" using the parabolic potential approximation
*/
double Grid :: calculateFrequency(Site::type from, Site::type to, Direction::type dir)
{
	double diffusionBarrier, prefactor;
	if(Direction::isVertical(dir))
	{
		diffusionBarrier = mBarrierVertical;
		prefactor = mFrequencyVertical;
	}
	else
	{
		diffusionBarrier = mBarrierHorizontal;
		prefactor = mFrequencyHorizontal;
	}
	double deltaEnergy = mEnergies[from] - mEnergies[to];
	double activationBarrier;
	
	if(deltaEnergy < 4 * diffusionBarrier)
	{
		activationBarrier = (4 * diffusionBarrier - deltaEnergy) * (4 * diffusionBarrier - deltaEnergy) / (16 * diffusionBarrier);
	}
	else
	{
		activationBarrier = 0.0;
	}
	Trace("Activation barrier from " << Site::getName(from) << " to " << Site::getName(to) << " direction " << Direction::getName(dir) << " = " << activationBarrier);
	return prefactor * exp(-activationBarrier / (mBoltzmannConst * mTemperature));
}

/*
	Counts the number of kinks.
*/
uint32_t Grid :: getKinkCount()
{
	uint32_t count = 0;
	
	for(auto i : mObjects)
	{
		if(i.second->getType() == Object::AlAtom)
		{
			GridAtom* atom = (GridAtom*)i.second;
			switch(atom->getSite())
			{
				case Site::Kink:
				case Site::DimerKink:
					count++;
					break;
				
				case Site::DoubleKink:
				case Site::DimerDoubleKink:
					count += 2;
					break;
			}
		}
	}

	return count / 2; //each kink was counted twice
}

void Grid :: removeAllEvents(GridAtom* atom)
{
	for(Direction::type dir = 0; dir < Direction::Count; ++dir)
	{
		removeEvent(atom, dir);
	}
}

void Grid :: removeEvent(GridAtom* atom, Direction::type dir)
{
	Site::type to = atom->getSite(dir);
	if(to != Site::Count)
	{
		Site::type from = atom->getSite();
		Trace("Removing hop event " << Direction::getName(dir) << " of atom " << atom << " from site " << Site::getName(from));
		auto hopRate = getRate(from, to, dir);
		mStateHash.togglePosition(hopRate->getHash());
		hopRate->removeAtom(atom);
		mStateHash.togglePosition(hopRate->getHash());
	}
}

void Grid :: addAllEvents(GridAtom* atom)
{
	
	for(Direction::type dir = 0; dir < Direction::Count; ++dir)
	{
		addEvent(atom, dir);
	}
}

void Grid :: addEvent(GridAtom* atom, Direction::type dir)
{
	Site::type to = atom->getSite(dir);
	if(to != Site::Count)
	{
		Site::type from = atom->getSite();
		Trace("Adding hop event " << Direction::getName(dir) << " of atom " << atom << " from site " << Site::getName(from));
		auto hopRate = getRate(from, to, dir);
		mStateHash.togglePosition(hopRate->getHash());
		hopRate->addAtom(atom);
		mStateHash.togglePosition(hopRate->getHash());
	}
}

GridAtom* Grid :: createAtom()
{
	if(mAtoms.size() == mAtoms.capacity())
	{
		throw Exception("Exceeded the reserved capacity of " + Utils::convertToString(mAtoms.capacity()) + " atoms");
	}
	mAtoms.push_back(AlAtom());
	return &mAtoms.back();
}

ostream& operator<<(ostream& os, Grid& g)
{
 	ostringstream builder;
 	uint64_t hash = UINT64_C(0);
 	
 	builder << "Size " << g.mRows << " " << g.mCols << endl;
	
	for(auto i : g.mObjects)
	{
		auto type = i.second->getType();
		uint32_t x = i.first.getx();
		uint32_t y = i.first.gety();
		uint64_t index = ((uint64_t)(y * g.mCols + x) << 8) | type;
		builder << type << " " << x << " " << y << endl;
		hash ^= Utils::getPrngKey(index);
	}
	string builderString = builder.str();
	os << "Checksum " << hex << setfill('0') << setw(16) << hash << endl << builderString;
	
	return os;
}

 		//     |D|
		//   |L|*|L| <-- even row
		//   |K|S|K|
		
		//   |K|S|K|
		//   |L|*|L| 
		//     |D|  <-- even row
Site::type Grid :: getSiteType(GridField& gf, GridObject* self)
{
	Trace("Getting site type of (" << gf.getx() << "," << gf.gety() << ") for object [" << self << "]");
	auto top = getObject(getTop(gf));
	auto bottom = getObject(getBottom(gf));
	auto left = getObject(getLeft(gf));
	auto right = getObject(getRight(gf));
	
	bool single, lateral, dimer, kinkLeft, kinkRight;
	GridObject *extraLeft, *extraRight;
	
	if(gf.isEvenRow())
 	{
		extraLeft = getObject(getBottomLeft(gf));
		extraRight = getObject(getBottomRight(gf));
		dimer = top != &GridObject::None && top != self;
		single = bottom != &GridObject::None && bottom != self;
	}
	else
	{
		extraLeft = getObject(getTopLeft(gf));
		extraRight = getObject(getTopRight(gf));
		dimer = bottom != &GridObject::None && bottom != self;
		single = top != &GridObject::None && top != self;
	}
	
	lateral = (left != &GridObject::None && left != self) || (right != &GridObject::None && right != self);
	kinkLeft = extraLeft != &GridObject::None;
	kinkRight = extraRight != &GridObject::None;
	
	if(dimer)
	{
		if(single)
		{
			return Site::DimerSingle;
		}
		else if(kinkLeft && kinkRight)
		{
			return Site::DimerDoubleKink;
		}
		else if(kinkLeft || kinkRight)
		{
			return Site::DimerKink;
		}
		return Site::Dimer;
	}
	else if(lateral)
	{
		return Site::Lateral;
	}
	else if(single)
	{
		return Site::Single;
	}
	else if(kinkLeft && kinkRight)
	{
		return Site::DoubleKink;
	}
	else if(kinkLeft || kinkRight)
	{
		return Site::Kink;
	}
	return Site::Free;
} 

inline void Grid :: setAtomSite(GridAtom* atom, Direction::type dir)
{
	auto gf = atom->getPosition();
	auto ngf = getField(gf, dir);	
	Site::type site;
	if(!isOccupied(ngf))
	{
		site = getSiteType(ngf, atom);
	}
	else
	{
		site = Site::Count;
	}
	Trace("Field in direction " << Direction::getName(dir) << ": (" << ngf.getx() << "," << ngf.gety() << ") = " << Site::getName(site));
	atom->setSite(dir, site);
}

void Grid :: refreshEvents(GridAtom* atom)
{
	Trace("Refreshing events for atom [" << atom << "]");
	removeAllEvents(atom);
	atom->setSite(getSiteType(atom->getPosition(), atom));
	for(Direction::type dir = 0; dir < Direction::Count; ++dir)
	{
		setAtomSite(atom, dir);
	}
	addAllEvents(atom);
}

void Grid :: refreshAtomEvents(GridObject * go)
{
	if(go->getType() != Object::None)
	{
		auto atom = (GridAtom*)go;
		Trace("Refreshing all events for atom *[" << atom << "]");
		refreshEvents(atom);
	}
}

void Grid :: refreshAtomEvent(GridObject * go, Direction::type dir)
{
	if(go->getType() != Object::None)
	{
		auto atom = (GridAtom*)go;
		Trace("Refreshing event " << Direction::getName(dir) << " for atom *[" << atom << "]");
		removeEvent(atom, dir);
		setAtomSite(atom, dir);
		addEvent(atom, dir);
	}
}

//affected positions:
		//       |0|
 		//     |8|1|9|
		//   |4|5|*|6|7| <-- even row
		//   |c|a|2|b|e|
		//     |d|3|f|
		
//affected positions:
		//     |d|0|f|
 		//   |c|8|1|9|e|
		//   |4|5|*|6|7| 
		//     |a|2|b|  <-- even row
		//       |3|
		
void Grid :: updateNeighbors(GridField& gf)
{
	Trace("Updating neighbors");
	
	refreshAtomEvent(getObject(getFarTop(gf)), Direction::Down);
	refreshAtomEvents(getObject(getTop(gf)));
	refreshAtomEvents(getObject(getBottom(gf)));
	refreshAtomEvent(getObject(getFarBottom(gf)), Direction::Up);
	refreshAtomEvent(getObject(getFarLeft(gf)), Direction::Right);
	refreshAtomEvents(getObject(getLeft(gf)));
	refreshAtomEvents(getObject(getRight(gf)));
	refreshAtomEvent(getObject(getFarRight(gf)), Direction::Left);
	
	auto x8 = getObject(getTopLeft(gf));
	auto x9 = getObject(getTopRight(gf));
	auto xa = getObject(getBottomLeft(gf));
	auto xb = getObject(getBottomRight(gf));
	
	if(gf.isEvenRow())
 	{
 		refreshAtomEvent(x8, Direction::Down);
		refreshAtomEvent(x8, Direction::Right);
		refreshAtomEvent(x9, Direction::Down);
		refreshAtomEvent(x9, Direction::Left);
		refreshAtomEvents(xa);
		refreshAtomEvents(xb);
		refreshAtomEvent(getObject(getBottomFarLeft(gf)), Direction::Right);
		refreshAtomEvent(getObject(getFarBottomLeft(gf)), Direction::Up);
		refreshAtomEvent(getObject(getBottomFarRight(gf)), Direction::Left);
		refreshAtomEvent(getObject(getFarBottomRight(gf)), Direction::Up);
	}
	else
	{
		refreshAtomEvent(xa, Direction::Up);
		refreshAtomEvent(xa, Direction::Right);
		refreshAtomEvent(xb, Direction::Up);
		refreshAtomEvent(xb, Direction::Left);
		refreshAtomEvents(x8);
		refreshAtomEvents(x9);
		refreshAtomEvent(getObject(getTopFarLeft(gf)), Direction::Right);
		refreshAtomEvent(getObject(getFarTopLeft(gf)), Direction::Down);
		refreshAtomEvent(getObject(getTopFarRight(gf)), Direction::Left);
		refreshAtomEvent(getObject(getFarTopRight(gf)), Direction::Down);
	}
}

void Grid :: removeAtom(GridAtom* atom)
{
	GridField gf = atom->getPosition();
	Trace("GridAtom *[" << atom << "] Removing from position: " << gf.getx() << "," << gf.gety());
	mObjects.erase(gf);
	updateNeighbors(gf);
}


void Grid :: placeAtom(GridAtom* atom, GridField gf)
{
	Trace("GridAtom *[" << atom << "] Placing at position: " << gf.getx() << "," << gf.gety());
	atom->setPosition(gf);
	mObjects.insert(FieldObjectPair(gf, atom));
	refreshEvents(atom);
	updateNeighbors(gf);
}

#ifdef GRIDCHECK
void Grid :: checkConsistency()
{
 	for(auto i = UINT32_C(0); i < mAtoms.size(); ++i)
 	{
 		checkNeighbors(&mAtoms[i]);
 	}
 	Trace("[GRID] Consistency check OK");
}

void Grid :: checkNeighbors(GridAtom* atom)
{
 	GridField gf = atom->getPosition();
 	auto site = atom->getSite();
 	auto realSite = getSiteType(gf, atom);
 	
 	if(site != realSite)
 	{
 		throw Exception("Atom *[" + Utils::convertToString(atom) + "] in position " +
				Utils::convertToString(gf.getx()) + "," + Utils::convertToString(gf.gety()) + ". Real site: " + 
 	 			Utils::convertToString(realSite) + ", Reported: " + Utils::convertToString(site));
 	}
 	
 	for(Direction::type dir = 0; dir < Direction::Count; ++dir)
 	{
		auto to = atom->getSite(dir);
		auto ngf = getField(gf, dir);
		if(to == Site::Count)
		{
			if(!isOccupied(ngf))
			{
				throw Exception("Atom *[" + Utils::convertToString(atom) + "] in position " +
					Utils::convertToString(gf.getx()) + "," + Utils::convertToString(gf.gety()) + ". Site in direction " +
					Direction::getName(dir) + " is not blocked");
			}
		}
		else
		{
			auto realTo = getSiteType(ngf, atom);
			if(to != realTo)
			{
				throw Exception("Atom *[" + Utils::convertToString(atom) + "] in position " +
					Utils::convertToString(gf.getx()) + "," + Utils::convertToString(gf.gety()) + ". Real site in direction " +
					Direction::getName(dir) + ": " + Site::getName(realTo) + ", Reported: " + Site::getName(to));
			}
			Trace("Checking rate for atom [" << atom << "] in position (" << gf.getx() << "," << gf.gety() << ")");
			auto rate = getRate(site, to, dir);
			if(atom != rate->getAtom(atom->getIndex(dir)))
			{
				throw Exception("Atom *[" + Utils::convertToString(atom) + "] in position " +
					Utils::convertToString(gf.getx()) + "," + Utils::convertToString(gf.gety()) + ". Not found in hop event " +
					Direction::getName(dir) + " from " + Site::getName(site) + " to " + Site::getName(to));
			}
		}	
 	}
 
}
#endif
